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Abstract 

The extinction of the contact process for epidemics in lattice models with quenched disorder 
is analysed in the limit of small density of infected sites. It is shown that the problem in such a 
regime can be mapped to the quantum-mechanical one characterized by the Anderson Hamiltonian 
for an electron in a random lattice. It is demonstrated both analytically (self-consistent mean- 
field) and numerically (by direct diagonalization of the Hamiltonian and by means of cellular 
automata simulations) that disorder enhances the contact process given the mean values of random 
parameters are not influenced by disorder. 
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I. INTRODUCTION 



The spread of epidemics in complex networks such as biological populations and computer 
networks is of great current interest, both for practical applications and from a fundamental 
point of view [J 3, 3, fl Q, El 

This is one of the issues of the theory of non-equilibrium phase transitions 00 



the theory of complex networks ^, |^ . The problems of interest include the question about 
the existence of a critical regime separating invasive (active) and non-invasive (absorbing) 
states of the system and, if such a transition exists, how it depends on internal and external 
parameters and also what the universal features of the transition are (see e.g. Q]). 

In one of the simplest models of epidemics, all the nodes are divided into two classes: 
infectious (/) and and susceptible (S) [8J]. The epidemic spreads by a contact process accord- 
ing to which an infected node can transfer infection to another susceptible node with typical 
infection rate w and recover with typical recovery rate e becoming again susceptible (the 
SIS model). The system undergoes a phase transition with variation of the dimensionless 
parameter, 



and 



IS 7] c 



= w/e, from the absorbing (rj < rj c ) to active state (rj > r] c ). The critical value 
,[7|, with Z being the typical number of links per node (coordination number). 



Usually, the infection and recovery rates are assumed to be node independent. However, 
in real systems, the values of w and e can vary from node to node (quenched disorder), 
investigations of contact processes in systems with quenched disorder over recent years 



1 fl 3, y Q Q, Q Q Q Q Q 



have resulted in some rather intriguing findings. 



For example, it has been suggested that the disorder can change the universality class of 



the model 



18, 



□ 



However, the situation is far from being completely understood, and 



the aim of this paper is to investigate the influence of a general form of quenched disorder 
on the dynamics of the contact process in the absorbing state. Using a combination of a 
simple epidemiological model with methods from condensed matter physics, we show how 
disorder in the infection or recovery rates, influences the long-time dynamics (decay time) 
of epidemics in the absorbing state. This is of practical importance in determining the time 
to extinction of epidemics within this state. We also identify a lower bound for r\ c and show 
how the degree of disorder influences the magnitude of the extinction rate. 

We consider the dynamics of the contact process far in the absorbing state when the prob- 
lem can be mapped to the quantum-mechanical one described by the disordered Hamiltonian 
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of the Anderson-type (see e.g. 20]) and an approximate method (self-consistent mean-field) 
can be applied. The spectrum of the Hamiltonian under this approximation is then used in 
the analysis of the long-time dynamics of the system. The advantage of the approach is in 
the possibility of incorporating a general type of disorder in the analysis while the disad- 
vantage is due to the rather severe restriction of being in the absorbing state (dilute regime 
for concentration of infected nodes). Our main result is that the disorder slows down the 
long-time dynamics of the system given the mean values of the random values stay the same 
as in ordered systems. The approximate analytical results are supported by exact numerical 
analysis using a cellular automata approach. 

The paper is organized in the following manner. The formulation of the problem is given 
in Sec. |HI The solutions in the dilute regime both for ordered and disordered cases are 



presented in Sec. IIHI followed by discussion in Sec. IIVI The conclusions are made in Sec. El 



II. FORMULATION OF THE PROBLEM 



Consider a set of N nodes (sites) connected to each other by links (infection paths). Each 
node, i, can be in one of two states: infected (occupied by an "excitation" and characterized 
by occupation number, n« = 1) or not infected (empty with n« = 0). The occupation 
number rii changes from to 1 as a result of infection from an occupied node j occurring 
with infection rate Wji, and from rii = 1 to rii = due to natural recovery with rate £j. Any 
state of the system is characterized by the set of occupation numbers, {n} = {rii, . . . ,n N }. 
Bearing in mind the stochastic nature of the infection and recovery processes it is convenient 
to characterize the system by the state vector \P(t)), the components of which are the 
probabilities of finding the system in different states at time t, \P(t)) = |P{ n }(i)), where 
n = 1, . . . , 2 N runs over all the possible states of the system. The time evolution of the state 
vector is governed by the master equation describing the conserved probability flow [3], 

8t\P(t)) = £\P(t)) , (1) 

where £ stands for the non-Hermitian Liouville operator, the non-zero elements of which 
describe the transitions between the states with different numbers of occupied nodes. 

It is convenient to make a linear transformation of the state coordinates (change of basis) 
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from P ni ,n 2 ,...,n N to the n-site probabilities, 

= ^ ] Pn 1 ,...,n k ,...,n i) ...,n N {t) > Pijif) = ^ ] Pn 1 ,...,n k ,...,n i ,...,n j ,...,n N if) > etc. , (2) 

where Pj(£) is the probability of finding node z in an occupied (n* = 1) state independent 
of the occupation of all other nodes. This allows the master equation ((TJ) to be recast in the 
following form: 

$P,(f) = -etPiit) +Y^w j i(P j (t) -Pji(t)) , (3) 

where Pj(t) — Pji(t) is the probability of finding the system with the occupied j-th node 
and the unoccupied i-th node, independent of the state of all other nodes. The single-site 
probability Pi(t) in Eq. Q is coupled with the double-site probabilities P%j{t). A similar 
probability-balance equation for Pij(t) contains the three-site probabilities and so on. This 
makes the set of simultaneous equations to be coupled and thus be non-trivial for analysis. 

The lowest level of approximations in decoupling schemes involves a complete ignorance of 
the double-site occupations, Py, in comparison with other terms in the master equation 
which is possible if 

P^^P^t) or Py(t) < —Pi(t) , (4) 

Wji 

for each pair of communicating sites % — j so that the master equation under these approxi- 
mations transforms to the following form, 

d t Pi(t) = -EiPi(t) + w jiPj(t) , ( 5 ) 

j 

and is hereafter referred to as the approximate master equation. The inequalities given by 
Eq. (0J) are valid if the typical recovery rate is much greater than the typical infection rate, 
e ^> w, and an epidemic dies out very quickly over the typical time, e^ 1 . In such a regime, 
the single-site probabilities are small for the majority of sites practically for all times, i.e. 
Pi -C 1, and this regime can be called a dilute regime for the concentration of infected sites. 

Therefore, making approximations given by Eq. (jJJ) we focus on the dynamics of the 
system far in the absorbing state ( rj rj c ), i.e. where an epidemic will certainly become 
extinct. Bearing in mind that the terms oc Py (entering Eq. (J3J) with a minus sign) reduce the 
infection rate due to a possible simultaneous occupation of both communicating sites % and 
j (the transmission of infection cannot occur between two sites if both of them are already 
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infected) we might expect that the solution of approximate rate equation © exhibits the 
enhancement of an epidemic in the dilute regime. In fact, the approximate master Eq. (JSJ 
on its own describes the spread of an epidemic in the system of N nodes with multiple 
reinfection of infected nodes where each node can be multiply occupied by the "excitations" 
(infection) and Pj has the meaning of an occupation number which can be larger than one. 

The approximate master equation (Eq. (jSJ)) similarly to the exact one can also have 
solutions which behave differently as t — > oo depending on the typical value of the parameter 
7]. In fact, there exists a critical value r)* for the approximate master equation which separates 
the absorbing and active states and this critical value rf c is smaller than the critical value r) c 
for the exact master equation due to the nature of the approximations made. This allows 
the lower bound estimate, rj*, for r] c to be found by solving the approximate problem for a 
quite general type of disorder. 

The state of the system in the dilute regime can be defined in the subspace of the singly- 
occupied sites spanned by the orthonormal site basis and is characterized by the state 
vector \P) = \Pi, P 2 , . . . , Pn) with components Pi(t) = (i\P(t)). The master equation (jHJ) 
can be recast as 

d t \P(t))=H\P(t)) , (6) 



where stands for the Liouville operator which now (under assumption of symmetric infec- 
tion rates, wa = Wji) is Hermitian and can be associated with the Anderson-like Hamiltonian 
(see e.g. \$), 

N 

in which the recovery rates, play the role of the on-site energies and the infection rates, 
Wij, can be associated with the transfer (hopping) integrals. Both of these values are random 
and this makes the further analysis non-trivial even in the dilute regime. The topology of the 
underlying network of sites, in principle, can be arbitrary but the simplest choice is a regular 
/^-dimensional lattice with nearest-neighbour interactions only. Furthermore, for simplicity, 
we consider a square lattice and thus the second sum in Eq. (J7J) runs for each site over 
Z = 4 its nearest neighbours only. It is worth mentioning that if the on-site matrix elements 
satisfy the sum rule, Si = ^ZjWij, then the number of occupied sites is conserved and the 
problem is equivalent to the random walk problem on a lattice with random transition 
rates j^jj, or to the scalar vibrational problem for a lattice with force-constant disorder 
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Notice also that decoupling of the single-site probabilities can also be made in the non-dilute 
regime by assuming that Pij = PiPj, i.e. ignoring possible correlations in occupation of 
the communicating sites. This brings a non-linearity to the problem which can be treated 
within the mean-field approach for ideal lattices []|. 



III. SOLUTION 



The formal solution of the problem given by Eq. (jSJ) is straightforward, 

\P(t)) = e**|P(0)> = ^e^V'lPtO))^) , (8) 

3 

where \eP) = |e{, . . . ,e° N ) and Xj are the eigenvectors and eigenvalues of the Hamiltonian, 
respectively, 7i|e J ) = A 3 -|e J '). Equivalently, this solution can be written via the Laplace 
transform of the state vector, \P(\)) = J °° \P(t))e~ xt dt, 

|P(A)> = (A - n)' 1 |P(0)> = g\p(o)) = £(A - A.rVi^o))^') , (9) 

3 

where Q = (A — 7^) _1 is the resolvent operator. 

We are interested in the time evolution of the total number of infected sites, 

J (*) = jf (E F '(*; , "o) ) = jf (flQFte*))) . ( 10 ) 

\ iio I \ Ho I 

averaged over different realizations of disorder (angular brackets) and/or over initial con- 
ditions (for concreteness, a single site i is infected at t — 0, i.e. Pi(0;i ) = 5u ) and its 
Laplace transform, 

J ( A ) = ^E^o(A)) • (11) 

Ho 

The other quantity of interest is the mean-squared displacement of the epidemic, 

i JL _ i JL 

(^))^E^< P >^»)> and ( R ^)) = jjlLKMioW) , (12) 

ii Ho 

where H ioi is the vector connecting site i with site i. 

As follows from Eq. (jHJ the dynamics of the system in the dilute regime are defined by 
the eigenspectrum of the Hamiltonian. The set of characteristic times (inverse eigenvalues) 
controls the evolution of the system with different eigenvalues being important for different 



time scales. The long-time dynamics of the system are defined by the maximum eigenvalue of 
the Hamiltonian, A max (A max < in the dilute regime) and our aim is to find an estimate for 
Amax and how it depends on the degree of disorder. The maximum eigenvalue depends on all 
the recovery and infection rates, A max (£j, and this obviously complicates its analytical 
evaluation. The exact analytical solution of the problem in the general case is not currently 
mown and numerous approximate analytical (see e.g. 0, Q) an d numerical (see e.g. 



24j,|25j) methods have been developed for evaluation of the spectrum of disordered Hamil- 



tonians. Below, we use one of the well-developed self-consistent mean-field approaches (th e 
homomorphic cluster approximation within the coherent potential approximation 2fil u\) 
to find the estimates for A max in the case of off-diagonal disorder and compare these results 
with the exact numerical calculations both for the Hamiltonian used in the approximate ap- 
proach and for the original problem (cellular automata (CA) calculations). Before analysing 
the disordered system, we start, however, with a trivial case of an ideal crystalline lattice in 
order to illuminate our approach. 



A. Ideal lattice 



In the case of an ideal crystalline lattice, the probability distribution functions are 5- 
functions, p e (£j) = S(si — s ) and p w {wij) = 8(wij — W), and the translationally invariant 
solutions of the eigenproblem are well known (see e.g.[28j). so that the eigenvectors are the 
Bloch's waves characterized by the wavevector k, 

|e k ) = N-V^e^li) , (13) 

i 

with Rj being the position vector of site % and the eigenvalues are 

A(k) = -e + w S k , (14) 

where = Ylj e _lkRi:7 ' (the sum is taken over j running over the nearest neighbours to 
arbitrary site i) is the structure factor and the wavevector k lies in the first Brillouin zone 
of the reciprocal space so that A max = Ak=o = —£o + Zwq. 

The Laplace-transform of the total number of infected states is then (see Eq. (fTTjl ). 
/(A) = (A — Ak=o) _1 ) and thus 

I(t) = e Ak =°* = e { - £o+WoZ)t . (15) 
7 



This means that in the dilute regime, r\ = Wq/eq <C 1, the exponent in Eq. ()15|) is negative 
and the total number of infected states decays exponentially with time. 

Therefore, for an ideal crystalline lattice the critical value of parameter rj = Wq/sq ob- 
tained from equation, A max = — 6q + Zwq = 0, for the system described by the approximate 
master equation is 

^ryst = , (16) 

which gives ^* ryst = 0.25 for a square lattice with nearest-neighbour interactions only. This 
estimate is equivalent to the standard (not self-consistent) mean- field estimate, and, as 
expected, is less than the true critical value, ^ crys t — 0.4122 7[. 

The Laplace transform of the mean-squared displacement (R 2 (X)) for ideal crystal in 
the dilute regime is given by the following expression (see Eq. (fT2"j) ). (R 2 (X)) = w a 2 Z(X — 
Ak=o) -2 , with a being the nearest neighbour distance, and thus 

(R 2 {t)) = w a 2 Zte Xk = ot = w a 2 Zte { - £0+WoZ)t . (17) 

It follows from Eq. (|17|) that the mean squared displacement increases exponentially in the 
active state when 77 > ^ryst (i- e - Vax = Ak=o > 0) and exponentially decays with time in 
absorbing state for r\ < r]* ryst (i.e. A max < 0). 



B. Disordered lattice 



The problems becomes much harder for a disordered lattice characterized by random 
infection and recovery rates. In order to find the time-dependence of the number of infected 
sites and the mean-squared displacement for the contact process we need to evaluate the 
configurationally averaged resolvent operator, (Q) (see Eqs. (fnjl- (fT2"j) ). This can be done ap- 
proximately for lattice models with certain types of disorder, namely, with diagonal disorder 
(disorder in the recovery rates, £j), off-diagonal disorder (disorder in transfer rates, 
and for binary s yste ms with substitutional disorder (two species of sites randomly occupy 
the lattice sites |2flj). One of the successful approximate analytical approaches is the self- 
consistent mean-field approach (coherent potential approximation, CPA) which allows the 
main spectral features of the disordered Hamiltonian and its eigenfunctions to be modelled 
(see e.g. 



S 



The main idea of the CPA is in replacement of the disordered lattice by the ideal crys- 
talline one which is characterized by the effective complex parameters (complex fields), e.g. 
by the effective recovery, e(X) = e'(X) + i§"(X), and transmission, w(X) = w'(X) + iw"(X), 
rates which depend on the eigenvalues, A, of the Hamiltonian and should be found self- 
consistently. The self-consistency equation follows from the requirement that a single defect 
placed in the effective crystal does not scatter the effective crystalline eigenfunctions if av- 
eraged over disorder. 

In what follows, for concreteness, we consider the case of off-diagonal disorder, when all 
the recovery rates are the same, p{ei) = 5{ei — eq), while the transfer rates are taken from a 
uniform (box) distribution, 

( (2A)- 1 if w - A < Wij < w + A 
P{Wij) = < , (18) 

I otherwise 

where A is the half- width of the distribution, < A < w , and the mean value Wij coincides 
with the crystalline one, w~ij = Wq. The particular form of the distribution (fTSj) is not 
important for the method discussed below. The conclusions are also applicable to any well 
behaved distribution given the mean value coincides with the value for the ordered system. 
The disordered Hamiltonian J7J) can be conveniently rewritten in the bond representation 



H = J2 (-Z-^iil - Z- X Ej\f)ij\ +^|i)01 +w ji \j){i\) , (19) 

where the summation is taken over all bonds (ij) in the system. Such a form of the Hamil- 
tonian allows the single non-correlated scatters (bonds) to be introduced in the absence of 

nn 

the on-site disorder (the homomorphic cluster approximation [2a, |2j| ) . The next step is to 
replace the above Hamiltonian with the effective non-Hermitian one, 

H = J2 {-Z- X E\i)(i\ - Z-H\j)(j\+w\i)(j\+w\j)(i\) , (20) 
(y) 

where the effective fields e and w are found from the following two self-consistency equations 
(see Appendix [XJ) [3^ . 

Z- 1 (e-e )±(w ij -w) \_ Q 
1 - {Q u ± %){Z-\e - e ) ± (wij -w))j 

The averaging in Eqs. (|2~Tj) is performed over random values of transition rates distributed 
according to the probability distribution given by Eq. (fTSj) . The effective resolvent (Green's 
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FIG. 1: (Color online) The spectrum of the effective (dashed line) and true (solid line) Hamiltonian 
(density of states) defined on the square lattice with nearest-neighbour interactions (Z = 4) in the 
system with transmission rates uniformly distributed around the mean value wq = 1 with half- 
width A = 1.0 and en = 30. The exact spectrum was obtained numerically using the kernel 



polynomial method 3J| for a model of N = 2000 x 2000 sites. The spectrum of the crystalline 
Hamiltonian with all nearest-neighbour interactions (A = 0) is shown by the dot-dashed line. The 
inset magnifies the spectrum around the top of the band. 

function) elements Gu and Gij can be expressed via the ideal crystalline resolvent elements, 
QH yst , of complex argument (see Appendix IXj) . 

1 



Gu(X) 



^O^cryst fWQ 



w 



w 



Zw 



(X + e)Gu(X)-l 



(22) 



which are well-known for the square lattice (see e.g. l 28J). 

The self-consistency equations fl21j) can be solved numerically and thus both effective 
fields can be found. Once the complex effective fields, e(A) and w(X), are known then the 
spectrum of the effective Hamiltonian can be found (see Fig. ^) enabling the dynamics of 
the system in the dilute regime within the self-consistent mean-field approach to be studied. 

It can be shown that the total number of infected states and the mean-squared displace- 
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ment in the CPA approximation obey the following equations (see Appendix iBj) . 

i fA max (A) Xt -r /-A max (A) At 

I(t) = Im / = dA = Im / -| — dA , (23) 

A min( A) A-A(A,k = 0) vr A min(A) X + e(X)-Zw(\) 

and 



1 /•A m ax(A) 

(i? 2 (t)) = Im / e xt 

^ Jx^JA) 



K </A min (A) 



VffiA.k) 
L(A-A(A,k))2j k=0 



dA 



1 /"^max(A) Za 2 w(A)e A * 
Im / ^ ? dA , 24 



7T 



with the effective dispersion law, 

A(A, k) = — e(A) + w(X)S-k ■ (25) 

The integration in Eqs. and (J23)l is performed over the band(s) of eigenvalues, A min (A) < 
A < A max (A), where the imaginary parts of the effective fields are finite for A > 0. 

As follows from Eqs. ()23|) and tfl^ the long-time dynamics of the system both for the 
number of infected nodes I(t) and for the mean-squared displacement (R 2 (t)) are defined 
by the largest eigenvalue. The upper band edge, A max (A), can be found within the CPA 
from the self-consistency Eqs. (J2*T|) by solving them for A > A max (A) where both effective 
fields are real, i.e. F±(w',e',\) = with F± standing for the left hand-side of Eqs. (|2~Tj). 
The analysis of the dependencies of the effective fields on A shows that the upper band edge 
corresponds to the branching point at which the following equation holds (see Appendix EJ) , 

~dF + dF_ dF_ OF. 



dw' de' dw' de' 



= . (26) 

- A max (A) 



The solution of Eq. (j2Ej) simultaneously with the self-consistency equations (|2T|) allows 
the position of the upper band edge, A max (A), to be found. The results of such an analysis 
are shown in Fig. |21 (see the solid line) for a particular choice of parameter t] ?7* ryst ~ 1 . 
When all the transfer rates are the same (A = 0), the maximum value coincides with the 
crystalline upper band edge, A max /w = —r]^ 1 + Z. When the disorder is introduced to the 
system the upper band edge shifts to larger values of A and thus the long-time dynamics 
slow down. The value of the shift increases with increasing degree of disorder characterized 
by the value of A (see Fig. |2J). This is a general effect that is independent of the type 
and form (probability distribution functions) of disorder given the mean values of random 
parameters are the same as for the ordered system. Indeed, any disorder brought to the 

11 




FIG. 2: (Color online) The dependence of the maximum eigenvalue, A max , evaluated using a mean- 
field approach (solid curve, labelled CPA), and A^ ax calculated by direct diagonalization (DD, 
circles for 4000x4000 and triangles for 2000x2000 lattices; the error bars represent the standard 
deviations of the distribution of the maximum eigenvalues), on the degree of disorder characterized 
by the half-width of the box distribution A for eo = 30 and wq = 1. The squares (labelled 
CA) represent the long-time decay rates, Aca> obtained by the CA simulations, each data point 
corresponding to approximately 5 x 10 10 runs, on a 5 x 5 lattice. The inset shows a version of the 
CA data scaled vertically to clarify the trend. 

system is equivalent to introducing additional interactions between the ordered eigenstates 
which unavoidably result in the level-repelling effect for the bare (crystalline) eigenstates 
32^ that is to the broadening of the spectrum. Therefore the disorder- induced slowing down 
of the dynamics of the contact process in the dilute regime is a general effect (if disorder 
does not influence the mean values of random variables) . 

It is known that the spectrum of the Hamiltonian obtained within the self-consistent 
mean-field approach usually reproduces very well the main features of the spectrum of the 
Hamiltonian (see Fig. excluding some special points like singularities (e.g. the mid- 
band singularity) and band edges (see the inset in Fig. |TJ). For the eigenstates around the 
band edges, the fluctuations of random parameters are essential and they lead to strong 
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localization of the eigenstates around the band edges. In fact, the true (non mean-field) 
density of states shows exponentially decaying band tails instead of sharp band edges typical 
for the mean- field crystals (see the inset in Fig. HJ. The mean-field approach does not take 
into account such fluctuations and thus the mean-field value of A max is only the (low-bound) 
estimate of the true maximum eigenvalue of the Hamiltonian, A^ ax . The values of Aj^ ax 
are random values depending on a particular realization of disorder. We have calculated 
numerically (using the Lanczos method) the distributions of the maximum eigenvalues for 
system of different sizes (up to N = 4000 x 4000 sites) and different values of A. The results 
for (A^ iax (A)) (averaged over 500 disorder realizations) are shown in Fig. El The averaged 
maximum eigenvalue depends on N (see Fig. |2J) apparently approaching some limiting value 
which is certainly less than the band-edge, — Eq + Zw, for crystal with all transfer rates equal 
to w = wo + A. 



IV. DISCUSSION 



The main finding from our analysis is that the disorder slows down the dynamics of the 
contact process in the dilute regime when the system is far in the absorbing state given 
that the mean values of the random parameters are the same as for the ordered system. 
This conclusion is supported by the exact solution for the contact process using the cellular 
automata (CA) simulations which were performed using a continuous-time algorithm similar 
to the n-fold way (see e.g. 33]). The CA simulations were used for evaluation of the inverse 
decay time for I(t — > oo) (see squares in Fig. EJ) which can be compared with A max (solid 
line) and A^ ax (circles and triangles in Fig. EJ). 

First, we compare the approximate (CPA) long-time decay rate (magnitude of the max- 
imum eigenvalue) with the exact (CA) one for an ideal (ordered) lattice. The dilute regime 
approximation enhances the epidemic and thus results in smaller long-time decay rates for 
I(t) and (R 2 (t)). Indeed, as follows from Fig. Elfor A = 0, the exact value A C a (— —26.215 
for a particular typical choice of parameters, Eq = 30 and Wq = 1) is smaller than the 
approximate one A max = — e + Zw = —26. When disorder is incorporated in the system 
the approximate decay rate decreases (eigenvalue increases) with increasing disorder (see 
the solid line in Fig. EJ)- The exact value of the decay rate in the disordered system also 
decreases (Aca increases; see the inset in Fig. EJ) with increasing disorder thus confirming 
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the tendency found in the dilute approximation, although the increase in Aca is appreciably 
smaller than the increase in A max and/or A^, ax . 

The value of the decrease in the decay rate with increasing disorder, being proportional 
to the width of the tail, |A max (A = 0) - A max (A = w )\/w ~ |A max (A = 0) - A^ ax (A = 
w o)\/ w o ~ 0.1, is naturally small due to the assumptions made, i.e. the system is far in 
the absorbing state. However the sign of the effect is important and it cannot be predicted 
by the standard (not self-consistent) mean-field analysis for the contact processes [3] if the 
mean values in the disordered system coincide with those for the ordered one. The other 
comment concerns the dependence of the effect on the parameters of the system. The 
broadening of the spectrum of the Hamiltonian does not depend on the mean recovery rate, 
Eq. This means that if Eq increases then the absolute magnitude of the effect within the 
CA treatment, A C a(A = 0) - A C a(A = w ) tends to A max (A = 0) - A^ ax (A = w ). If the 
recovery rate decreases the system approaches the active state and approximations made 
for the dilute regime break down. Therefore the analysis performed cannot be considered 
as a good approximation around criticality. In fact, as it follows from the preliminary CA 
analysis, the value of Aca (A = 0) — Aca (A = wo) increases with decreasing e and can 
reach zero around criticality and even change the sign (to be discussed further elsewhere) 
i.e. the introduction of disorder into a crystalline system at criticality causes a transition to 
the absorbing state rather than further into the active phase. 

The last comment concerns a possible rough estimate of the critical parameter r] c for 
transition from the absorbing to active state in disordered system. This estimate can be 
found by solving the equation A max (A) = (or A^ ax (A) = 0), which gives, r\ c ~ ^ ryst (l + 
(A ma x(A = 0) — A max (A = w ))/w ). Of course, the quality of this estimate is the same as 
that for the crystal, i.e. r/£ ryst = Z~ l = 0.25 as compared to exact value r] c ~ 0.4122, and 
can serve only as a reliable low bound for the critical value. 

V. CONCLUSIONS 

We have presented the analysis of the contact process in the limit of low density of 
occupied (infected) sites (see Eq. Q), i.e. in the dilute regime for the infected sites when 
all the correlation effects in occupation probabilities can be ignored. This limit occurs, e.g. 
when the transfer rate is much smaller than the infection rate, w <C e. The system resides 
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in the absorbing state for such a range of parameters and its dynamics can be described by 
using the quantum-mechanical tight-binding Hamiltonian. The disorder, both in the transfer 
and recovery rates, can be incorporated into the formulation which can be reduced to the 
eigenproblem for the Anderson-like Hamiltonian (see Eq. (J7J)). The eigenproblem can be 
solved approximately analytically (self-consistent mean-field) and exactly numerically and 
thus the estimate of the decay rate, |A max |, for long-time dynamics can be found for different 
degrees of disorder (see Fig- EJ) • The approximate solution is supported by exact numerical 
analysis using the cellular automata approach. In particular, we conclude that any type 
of disorder which does not change the mean values of random parameters slows down the 
long-time dynamics of the contact process occurring in the far absorbing state. 



APPENDIX A: SELF-CONSISTENCY EQUATION WITHIN THE HOMOMOR- 
PHIC CLUSTER CPA 

In this Appendix, we derive the matrix self-consistency equation within the homomor- 
phic cluster CPA. Within the self-consistent mean-field approach the disordered lattice is 
replaced by an ideal lattice characterized by self-consistently found effective complex param- 
eters (fields), namely, by the effective recovery and transmission rates, e(X) and w(X). The 
effective Hamiltonian describing this effective lattice is given by Eq. (|2U|) . The effective fields 
are found within the single-defect approximation according to the following standard proce- 
dure 3,0- 

A single defect bond, (ij), taken from the random set of bonds characterizing 
the disordered lattice is placed in the effective medium. The single-defect Hamiltonian, Hi, 
is the sum of the ideal effective Hamiltonian and the perturbation due to the defect bond, 
n 1 = H + 5V, where 

SV = -Z-\e - e(X)) + + K - 0(A)) + , (Al) 

This defect bond influences (scatters) the eigenfuctions of the original effective Hamiltonian. 
In the CPA, the effective medium is tuned in such a manner that this scattering vanishes 
on average. In other words, the single-defect scattering operator T, introduced by equation 
T = 5V + 5VQT (where Q = (A — H)^ 1 is the effective resolvent) averaged over different 
realizations of defect bond taken from the same probability distribution as for disordered 
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medium should be zero, i.e. 

(f) = (5V(l-g5V)- 1 ) = , (A2) 

The above equation is the self-consistency matrix equation where the scattering matrix in 
the site basis is 

f _ 1 / 5e + (5w 2 - 5e 2 )Gjj 5w - (5w 2 - 5e 2 )G tJ \ 
~ \l-gSV\ y5w- (5w 2 - be 2 )Q 3l be + (6w 2 - 6e 2 )G n ) ' 

with 5e = —Z~ 1 (eq — e(A)) and 5w = Wij — w(X). Bearing in mind that Gu = Gjj and 
Gij = Gji the scattering matrix can be easily diagonalized by a similarity transformation to 
a new basis, so that 

8e—8w 







ji _ i.-(yn-yii)(oe-ow) (A4) 

n 8e+5w 1 

which straightforwardly results in Eq. (|21|) . 

The elements of the scattering matrix T depend on the diagonal and off-diagonal ma- 
trix elements of the effective resolvent, Gu and Gij, respectively. The effective Hamiltonian 
describes the ideal lattice characterized by complex parameters, e and w, and thus its eigen- 
fuctions are the Bloch's waves given by Eq. (|13|) with effective dispersion described by 
Eq. (|25jl . This allows the real-space matrix elements of the resolvent to be expressed via 
the reciprocal-space ones, and then via similar elements of the resolvent for ideal crystalline 
lattice, 



A" ^— ' A + e — wSk wN ^ Xw/w + ew/w — e — A cryst (k) 
w(X + e) 



^ ^->cryst 

W W 

■>cryst / 



^0 



(A5) 



with Gu ys (A) = A" 1 X/k(*^ — ^cryst(k)) 1 being the crystalline resolvent characterized by 
the crystalline dispersion, A cryst (k), given by Eq. (fPlj) and 

a (X) = 1 V e ' ikR ' J 1 W A + £ ~ 

A ' N ^ A - A(A, k) wZN ^ VA - A(A, k) 

= A±£I E _^ J_, (A6) 

wZ N ^ A — A(A, k) wZ ' y J 

where A^ 1 Xlk(^ — ^(^> = Ga(X). Eqs. (|A5|) - f|A6|) justify the expressions for the matrix 
elements of the effective resolvent given by Eq. 
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APPENDIX B: NUMBER OF INFECTED SITES WITHIN THE CPA 



The aim of this appendix is to derive Eq. (|23|) for the number of infected sites as a function 
of time. It is convenient to find the Laplace transform, /(A), of the function I(t) and then 
use the inverse transform, 



I(t) = — / e A */(A)dA , (Bl) 

to reveal Eq. ()23j) . The Laplace transform J(A) is given by Eq. (jllj) in which, within the 
CPA, the averaged resolvent matrix element should be replaced by the matrix element of 
the effective resolvent, 

1 JV 1 N 1 ikR 

, ( A) = _£^(A)=_£_£^J^. (B2, 

Bearing in mind the identity ^\e~ lkR "o = iV7) k) o, we obtain /(A) = (A — Ak =0 ) _1 . Substitu- 
tion of this expression in Eq. (jBl)) gives 



i /-(5+ioo \t 

I(t) = — / dA . (B3) 

Anl Js-ioo A — A k=0 

The effective dispersion, Ak, depends on the the effective fields w(X) and e(X) which are 
analytic functions of A everywhere except the finite interval on the real axis (branch cut), 
A G [A min , A max ], where the density of states of the effective Hamiltonian is finite (see e.g. 
j^J and references therein). The contour of integration in Eq. (1B3j) can be transformed into 
a closed one around the branch cut and thus, taking into account that the real part of the 
integrand is a continuous function through the branch cut but the imaginary part changes 
sign, Eq. (}2~3*j) follows from Eq. (jB3|) and Eq. ()24j) can be derived in a similar fashion. 



APPENDIX C: EQUATION FOR THE BAND EDGE 

The two self-consistency Eqs. (121 J) can be recast in the following form: 

5e — 5w 



-piwi^dwij = F_(f?,w,A) = (CI) 



o 



1 - (Gu - Qij)(Se - 8w) 

, £ t ™ piw^dwij = F+(e, w, A) = . (C2) 

'o 1 - {Gu + Qij){de + dw) 

For any value of A in the complex plane, these equations can be solved and two complex 
fields, e(A) and w(X), can be found. On the real axis for A, A > A max , the both fields are 
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also real, e(A) = s'(X) and w(X) = w'(X) with A max being the branching point. In this range 
of A, it is convenient to rewrite Eqs. ([C1I) - (1C2I) in the following form: 



e' — e-(w',X) 

, (C3) 

£ = £+[W , A) 

where e T are multivalued functions of w' for fixed A. If A > A max these contour lines usually 
cross at two points one of which corresponds to the physical solution. At the branching 
point, A = A max , these two solutions merge and the condition for this is 

dw' dw' ' ^ 
which together with Eqs. (jCl|) - (|C2|) or with Eqs. (jC3|) allows the location of the upper band, 
Amax; to be found. Eq. (jC4|) can be rewritten in the more elegant but equivalent form given 
by Eq. (|2*BI). Indeed, differentiation of Eqs. (|C1J) - (|C2J) with respect to w' gives 



de T dF T fdF T 



dw' dw' \ de' 
from which Eq. (J2*rjj) follows straightforwardly with the use of Eq. (|C4|) . 



(C5) 
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